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Abstract. The nature of transverse instabilities to dark solitons and dispersive shock waves for the (2+1)- 
dimensional defocusing nonlinear Schrodinger / Gross-Pi taevskii (NLS / GP) equation is considered. Special atten- 
tion is given to the small (shallow) amplitude regime, which limits to the Kadomtsev-Petviashvili (KP) equation. We 
study analytically and numerically the eigenvalues of the linearized NLS / GP equation. The dispersion relation for 
shallow solitons is obtained asymptotically beyond the KP limit. This yields 1) the maximum growth rate and as- 
sociated wavenumber of unstable perturbations; and 2) the separatrix between convective and absolute instabilities. 
The latter result is used to study the transition between convective and absolute instabilities of oblique dispersive 
shock waves (DSWs). Stationary and nonstationary oblique DSWs are constructed analytically and investigated nu- 
merically by direct simulations of the NLS / GP equation. The instability properties of oblique DSWs are found to 
be directly related to those of the dark soliton. It is found that stationary and nonstationary oblique DSWs have the 
same jump conditions in the shallow and hypersonic regimes. These results have application to controlling nonlinear 
waves in dispersive media. 

1. Introduction. The instability of one-dimensional structures to weak, long wave- 
length, transverse perturbations plays an important role in multi-dimensional nonlinear wave 
propagation. Examples include nonlinear optics [1], Bose-Einstein condensates (BECs) [2], 
and water waves [3, 4]. Early theoretical work on the transverse instability of solitons for 
the Kadomtsev-Petviashvili (KP) equation [5, 6] and the nonlinear Schrodinger (NLS) equa- 
tion [7, 8] focused on its existence and maximum growth rate, both properties of the real 
portion of the spectrum of unstable modes. Recent numerical simulations of NLS [9] and 
vector NLS [10] supersonic flow past an obstacle in two-dimensions reveal the excitation of 
apparently stable, oblique spatial dark solitons for certain flow parameters. The resolution of 
this inconsistency was explained in [11], where the instability was shown to be of convective 
type so that transverse perturbations are carried away by the flow parallel to the soliton plane, 
effectively stabilizing the soliton near the obstacle. The characterization of convective versus 
absolute instability requires knowledge of the spectrum for a range of wavenumbers in the 
complex plane [12, 13]. For NLS dark solitons, the criteria can be simplified and involve the 
imaginary (stable) portion of the spectrum [11]. 

One of the hallmarks of supersonic flow is the formation of shock waves. In classical, 
viscous fluids, shock dynamics can be well understood mathematically in the context of a dis- 
sipative regularization of conservation laws (cf. [ 14]). There are, however, a number of fluids 
with negligible dissipation whose dominant regularizing mechanism is dispersion (see the re- 
view [15]). Most notably, superfluidic BECs and optical waves in defocusing nonlinear media 
fall within this class of dispersive fluids. When a dispersive fluid flows at supersonic speed, it 
can form a dispersive shock wave (DSW) that possesses an expanding, oscillatory wavetrain 
with a large amplitude, soliton edge and small amplitude sound wave edge. DSWs appear 
as special, asymptotic solutions of nonlinear dispersive equations and have been observed in 
BEC [16, 17, 18] and nonlinear optics [19, 20]. Their theory is much less developed than 
their classical (dissipative) counterparts. In particular, there has been limited study of DSW 
stability. Recent works numerically observe transverse instabilities for NLS DSWs resulting 
from dark pulse propagation on a background in two spatial dimensions [21] and for oblique 
DSWs in supersonic flow past a corner [22]. In the former case, the transverse instability 
was mitigated by introducing nonlocal nonlinearity while in the latter case, the convective 
nature of the instability effectively stabilizes the oblique DSW in certain parameter regimes. 
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In contrast, oblique shock waves in multidimensional, classical gas dynamics are known to 
be linearly stable when the downstream flow is supersonic [23, 24, 25] (see also the review 
article [26] for more general results). 

The aim of this work is to clarify the role of absolute and convective instabilities as 
they relate to spatial dark solitons and apply this understanding to DSWs in multiple spatial 
dimensions. Analytical and computational challenges include: 

• The multi-dimensional nature of the flows. 

• The general criteria for absolute and convective instabilities requires detailed knowl- 
edge of the spectrum. 

• Long time integration and large spatial domains are required to properly resolve 
DSWs numerically. 

To address these challenges, we asymptotically determine the spectrum of transverse 
perturbations to shallow but finite amplitude NLS dark solitons beyond the Kadomtsev- 
Petviashvili (KP) limit. This enables determination of the maximum growth rate and associ- 
ated wavenumber of unstable perturbations. Using adjoint methods, we introduce a simple, 
accurate method for computating the spectrum and its derivatives numerically for arbitrary 
soliton amplitudes. Simplified criteria for the determination of the separatrix between ab- 
solute and convective instabilities are derived. The separatrix is determined in terms of the 
critical Mach number M CI as it relates to the soliton far field flow. Oblique dark solitons 
are convectively unstable when M > M CI and absolutely unstable otherwise. Using our 
asymptotic and numerical computations of the spectrum, we determine M cr , demonstrating 
that 1 < M cr ~ 1.4374 with M cr a monotonically increasing function of soliton amplitude. 

The oblique DSW trailing edge is well-approximated by an oblique dark soliton. In this 
study, we apply the soliton stability results to the oblique DSW trailing edge in the stationary 
and nonstationary cases. Stationary oblique DSWs result from the solution of a boundary 
value problem (supersonic corner flow) while the nonstationary case arises in the solution of 
a Riemann initial value problem. We find that oblique DSWs with supersonic downstream 
flows can be absolutely unstable in contrast to classical oblique shocks. We also show that 
stationary and nonstationary oblique DSWs have the same downstream flow properties in the 
shallow and hypersonic regimes. 

We consider the (2+l)-dimensional defocusing (repulsive) nonlinear Schrodinger/ 
Gross-Pitaevskii (NLS / GP) equation 

itk = -^(V»x* + ^ ra ) + MV , (x,y) GM 2 , t>0, (1.1) 
along with appropriate initial and/or boundary data. Equation (1.1) models matter waves 
in repulsive BECs and intense laser propagation in optically defocusing (i.e., with normal 
dispersion) media. In the variables 

i> = ^pe i *, (u,v)=V<t>, (1.2) 
Equation (1.1) can be recast in terms of the fluid-like variables p (density) and (u, v) (super- 
fluid velocity) 

pt + (pu) x + (pv)y = , (1.3a) 
u t + uu x + vuy + Px = \( Pxx ± Pvv , (1.3b) 



n + uc, + rr„ : = \ [ P** + Pw _ &+A | . (L3c) 



Note that eqs. (1.3) in the dispersionless regime (neglecting the right hand sides) correspond 



DARK SOLITONS, DSWS, AND TRANSVERSE INSTABILITIES 



3 



to the classical shallow water equations (Euler equations of gas dynamics with adiabatic con- 
stant 7 = 2) with the speed of sound y/p [27]. 

The outline of this paper is as follows. Section 2 discusses the spectrum of unstable trans- 
verse perturbations of dark solitons with asymptotic resolution of the maximum growth rate 
and associated wavenumber in the shallow regime. Using analytic properties of the spectrum, 
we recap the derivation of the general criteria for absolute and convective instabilities and for 
oblique solitons, we derive the simplified criteria in Sec. 3. The separatrix M CT is determined. 
We derive nonstationary oblique DSWs of arbitrary amplitude and stationary oblique DSWs 
in the shallow regime, showing the connection between their downstream flows in Sec. 4. The 
stationary case is compared with (2+l)-dimensional numerical simulations. Convective and 
absolute instability of oblique DSWs is described in terms of the separatrix for the trailing 
edge dark soliton. Our numerical methods are presented in Sec. 5. Finally, Sec. 6 contains 
a discussion of the results and the applicability of our methods to other nonlinear dispersive 
problems. 

2. Transverse instability of dark solitons. It is well-known that dark soliton solutions 
of (1.1) exhibit an instability to perturbations of sufficiently long wavelength in the transverse 
direction along the soliton plane [8] . The eigenvalue problem associated with linearizing (1.1) 
about the dark soliton leads to the dispersion relation for unstable perturbations. Beyond 
demonstrating the existence of an instability, knowledge of the dispersion relation for a range 
of wavenumbers yields important properties of the instability, such as the growth rate r max , 
the maximally unstable wavenumber fc max , and whether or not the instability is convective or 
absolute. 

An example numerical computation of the eigenvalues for the spectral problem in eq. (2.7)| 
is shown in Fig. 2.1. Since exact expressions are not known, asymptotic approaches lever- 
aging the shallow dark soliton, KP limit [6, 28] and others [29, 11] have been devised. In 
this section, we complement these results by determining the next order correction to the dis- 
persion relation for shallow dark solitons resulting in accurate approximation across a wider 
range of soliton amplitudes. We use this to determine r max and fc max asymptotically. These 
calculations are verified numerically. 




FIG. 2.1. The real (dashes) and imaginary (solid) parts of the discrete eigenvalue f!o(fc; f) of the linearized 
NLS equation (spectral problem (2.7)) as functions of k for v = 0.5. Delineated on the axes are: i) the cutoff 
wavenumber fc cu toff [Eq. (2.11)7, where the eigenvalue transitions from purely imaginary to real, ii) the maximal 
growth wavenumber and growth rate (fcmax, Tmax) [Eq. (2.14)7, ond Hi) the critical wavenumber and associated 
eigenvalue (fc C r, f2cr) [Eq. (3.8)7 corresponding to the transition between absolute/convective instabilities. 
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2.1. Dark Soliton. Up to spatio-temporal shifts and an overall phase, the most general 
line dark soliton solution of (1.1) is 



ip' s (x', y', t') = yfp{cos(f> + isin^tanh [a(sm(3x' — cos fly' — vt')]} (2.1) 

„2 , j2 



x exp < i 



cx' + dy' - [ ° ^ d +p)t' 



a = ^fp sin <f), v — c sin /3 — d cos (3 — yfp cos < 



■e [o,tt] , 



where p is the background density and the phase jump across the soliton 2<fi determines the 
depression amplitude as ^fp\ sin0|. The soliton is propagating at an angle f3 with respect 
to the (horizontal) x' axis, with horizontal and vertical flow velocities c and d, respectively. 
Interpreting this solution in the fluid context with density IV^I 2 and flow velocity V axgil>' s , 
the soliton is a localized density depression on a uniformly flowing background. The Mach 
number of the background flow is the total flow velocity divided by the speed of sound 



M = 



+ d 2 



(2.2) 



The soliton has the far field behavior 



sin f3x' — cos fly' — > ±oo , 



cx' + dy' — 



+ P)t' 



Thus, five parameters determine the soliton uniquely, i.e. p, <fi, fl,c,d. 

Using the invariances of Eq. (1.1) associated with rotation, Galilean transformation, scal- 
ing, and phase, we apply the coordinate transformation 



, z jX 2 —i \( c si n 3—d cos /3) -4= + {c cos B+d sin 8) -7= + 

ip s (x, y,t) = — e L VP vp 

Vp 



y 1 (^+d')t 



(2.3) 



ip' ( -^—(sinflx + cos fly) + —t, —=(— cosflx + sin fly) + —t, — ) 
\VP P \[P P p) 



leading to the one-parameter family of dark solitons 

V> s (£, y,t; v) = [in + ^tanh e~ lt , v 2 + n 2 = 1 , 
where v = | sin</>| £ (0, 1] and the frame moving with the soliton is 



(2.4) 



£ = x — Kt . 



The soliton amplitude is v. When v -C 1 the dark soliton is in the shallow amplitude regime. 
The soliton speed is n = — cos <j> = VT — v 2 . 

2.2. Linearized eigenvalue problem. To study the transverse instabilities of the dark 
soliton (2.4), we consider the ansatz forEq. (1.1) 



ip(€, v, *; v) = [ips{£, y, t)c lt + <p R (£, y, t) + , y, t)] 
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where ipR, ipj are the real and imaginary parts of a small perturbation. Linearizing (1.1) 
results in the system 



d_ 

at 



tp = Lip 



C = 



(2.5) 



2^fctanri(^£) 
-i/ 2 [2-3sech 2 (>£)] 



+ d vv 
kO^ — 



) - ^ 2 [2 + scch 2 K)] 
2vk tanh(y£) 



It is expedient to decompose the perturbation as 

Substituting (2.6) into (2.5) yields the linearized spectral problem 

J£f(£;fc) - -ifi(fc)f(£;fc), 

where 



A(£;fc) 
/a (6*0 



L = Lq + \k 2 



J 



and 



1 

-1 



-i% + ^ 2 [2-3sech 2 K)] 
nd^ + 2^Ktanh(i/^) 



—ndf 
-1% 



U>«-^2 



2i/k tanh(^£) 
sech 2 (i/f)] 



(2.6) 
(2.7) 
(2.8) 

(2.9) 
(2.10) 



For fc e R, L and i are self-adjoint with respect to the i 2 (M) inner product 

(g,h) = J R S T h*d(. 

For small k it was shown formally in [8] that: (;') a double eigenvalue f2(0) = bifurcates 
into two distinct branches with each in iR; (ii) there is another zero eigenvalue at the cutoff 
wavenumber 



fccutoff = - 2 + - V 2 + 1 , fi(fccutoff) = . (2.11) 

These calculations were made rigorous in [30] and can be summarized as follows. 

Theorem 2.1 (Rousset, Tzvetkov [30]). For k e (— fc cuto ff , fccutoff) \ {0}, the sys- 
tem (2.7) has exactly two purely imaginary eigenvalues which are simple and come in pairs 
=bOo(fe). Therefore, the dark soliton is unstable to sufficiently long wavelength transverse 
perturbations. Furthermore, for k € K, |fc| > fccutoff, the spectrum fi(fc) is real. 

For the study of convective/absolute instabilities, knowledge of the stable portion of the 
spectrum when |fc| > fc cu toff is required. Based on numerical and asymptotic computations, 
we conjecture the following. 

CONJECTURE 2.2. For \k\ > fccutoff, there exist exactly two real, simple eigenvalues 

±n (k). 

This conjecture is a natural extension of Thm. 2. 1 . See Appendix B for further details 
and comments. 

Without loss of generality, we choose f2o(fe) such that 3{f2o(fc)} > for < k < 
^cutoff and 5i{J7o(fc)} > for k > fccutoff- Thus, fio(fc) is the dispersion relation for trans- 
verse perturbations of the dark soliton (2.4). By suitable choice of a branch cut, the eigenvalue 
flo(k) can be analytically continued for k E C \ {0, ±fc cu toff} with and ±fc cuto ff square 
root branch points. We denote the growth rate as 

r(fc) = 3{ft (fc)} , (2-12) 



6 



M. A. HOEFER AND B. ILAN 



and the eigenfunction associated with Q (k) as 

fo(£;fc) = 



/o,i(6*0 



In Section 5 we discuss our numerical method for computing ilo(k) for k g C. To illus- 
trate the spectrum, Figure 2. 1 presents the dependence of the (real or imaginary) eigenvalue, 
D,o(k), on k. Figures 2.2 and 2.3 present the computed continuous and discrete spectra for 
particular wavenumbers < k < fc C utoff (^o € iM.) and k > fc C utoff (^o G R) as well as the 
associated localized eigenfunctions. Note that the eigenfunctions are neither symmetric nor 
anti-symmetric. 



0.1 



%{n} o 



-0.1 



(a) 




-22 



(c; 




22 



FIG. 2.2. (a) Numerical approximation of the continuous spectrum (•) and the two purely imaginary discrete 
eigenvalues ±f2o ~ ±0.022i (+) computed for the linearized system (2.7) with v = 0.5 and k = 0.2 < fc cu t ff ^ 
0.23 [Eq. (2.11)7. The real (solid) and imaginary (dashed) parts of the corresponding two component localized 
eigenfunction are shown in (h) /o,l(£) and (c) /o,2(£). 
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FIG. 2.3. Same as for Fig. 2.2 except k = 0.25 > fc cuto ff and ±S1q ~ ±0.0215. 



2.3. Asymptotic eigenvalue. It follows from (2. 1 1) that for shallow amplitude solitons, 
< v <C 1, the cutoff wavenumber is small, i.e., fc cuto ff ~ ^ y2 - ^ n Appendix A we prove: 

PROPOSITION 2.3. For shallow amplitude, < 1/ <C 1, and either k < fc C utoff or 
^cutoff < k ~ C(^ 2 ) "C 1, f/ie eigenvalue for (2.7) satisfies 

. 3 6V2-y/3fc - 3^ 2 



(2.13) 



KP 



NLS correction 

vv/zere the first (leading order) term is the dispersion relation for the KP equation and the 
second term is the 0(k 5 ^ 2 ) correction arising from the NLS equation. 
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Equation (2. 13) gives an asymptotic approximation to the eigenvalue for long wave per- 
turbations of shallow dark solitons. The dispersion relation for the KP equation is well known 
(cf. [6, 28]). The new 0(fc 5 / 2 ) correction term enables us to accomplish the following. 

• Implement an accurate, explicit calculation of the maximum growth rate and associ- 
ated wavenumberof unstable perturbations (Sec. 2.4). 

• Show that the separatrix between absolute and convective instabilities is supersonic 
(Sec. 3.3). 

• Validate the numerical computations of ilo(k), which are sensitive and computa- 
tionally demanding, especially in the shallow regime. 

2.4. Calculation of the maximum growth rate. The maximal growth wavenumber 
fcmax and the maximum growth rate r max are defined by 

fto(fc max ) = o, r max = 3{r> (fc m ax)} ■ (2- 14 ) 

Since D,o(k) is real for k > k cuto g it follows that & max < fc cu toff (see Fig. 2.1). Using 
Proposition 2.3 we find 
Corollary 2.4. 

* m "=i + H +0< " 8) ' <ll5) 

KP NLS correction 

+0(v 7 ). (2.16) 



3V3 9^3 

KP NLS correction 

The proof follows by expanding fc max and r max for small v and solving eq. (2.14) with the 
approximation (2. 13). 

A comparison of these results with numerical computations (discussed in Sec. 5) is 
shown in Fig. 2.4. The computations exhibit excellent agreement with the asymptotics as 
well as the expected scaling of the errors with v. 

3. Convective and absolute instabilities of dark solitons. We begin by reviewing the 
notions of absolute/convective instabilities and the general criteria for distinguishing between 
them. For more detailed discussions see [12, 13, 31, 32, 33]. Qualitatively, absolute and 
convective instabilities can be defined as follows (see illustration in Fig. 3.1). 

DEFINITION 3.1. A solution is said to be absolutely unstable if generic, small, localized 
perturbations grow arbitrarily large in time at each fixed point in space. A solution is said to 
be convectively unstable if small, localized perturbations grow arbitrarily large in time but 
decay to zero at any fixed point in space. 

It is important to note that Definition 3.1 depends implicitly on the reference frame as 
can be gleaned from Fig. 3. 1 where panel (b) is a rotation in the x-t plane of panel (a). Such 
a rotation implies that the observer in (b) is moving faster to the left than the observer in (a). 
Thus, if the observer "outruns" the growing perturbation, then the instability is convective. 
Equivalently, if the background flow speed is faster than the expanding, unstable perturbation, 
and after sufficient time passes the solution returns to is unperturbed state, the instability is 
convective. 

3.1. Review of the general criteria for distinguishing between instabilities. Abso- 
lute and convective instabilities can be distinguished analytically. Consider an initial value 
problem on the entire line, i.e. a (l+l)-dimensional linearized system on (x, t) £ Kx (0, oo). 
The usual approach for studying instabilities is to consider a small, spatially extended plane 
wave perturbation e l ( kx ~ ut ) G f some wavenumber k and corresponding frequency oj = il(k) 
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FIG. 2.4. Numerically computed maximum growth rate r max (a) and maximally unstable wavenumber fc ma x 
(c) as functions of v for dark solitons of the NLS equation (1.1). The KP limit and its first order correction are 
presented for comparison. Plots (b) and (d) are the corresponding differences between the highly-accurate computed 
values and asymptotic approximations (2.15)— (2.16), exhibiting the expected scaling with v. 



(a) (b) 




FIG. 3.1. Illustration of (a) absolutely and (b) convectively unstable waves. 



determined by a zero of the dispersion function T)(ui,k) = 0. The zero state is stable if and 
only if 3{f2(fc)} < for all zeros of the dispersion function. However, the evolution of a 
particular, localized perturbation involves a Fourier integral over all real wavenumbers so that 
treating a single wavenumber is insufficient to fully describe any instabilities observed (or not 
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observed) in a physical system [12]. 

The resolution calls for a different approach to instability analysis. Instead of a plane 
wave perturbation one assumes that the system is perturbed by a localized impulse, i.e. a 

Dirac delta function at position x and t = 0. In this case the solution is the Green's function 

r r pi(kx—ujt) 

G( X ,t) = JJ^y^dk, (3.1) 

where the Fourier integral is carried out over real wavenumbers and the Bromwich frequency 
contour lies above all zeros of T>(k,u). In connection with the plane-wave analysis, the 
system is unstable if and only if the solution grows without bound along some reference 
frame, i.e. there is a velocity V such that for fixed x, 

G(x — Vt, t) — H? oo unstable . 

However, when considering a particular reference frame, say £ = x — Vot for fixed Vq, if the 
solution grows without bound (resp. decays to zero) at a certain fixed point in space, x, then 
the system is absolutely (resp. convectively) unstable in this reference frame, i.e. 

• G(x — V)t, t) oo absolutely unstable. 

• G(x — Vot, t) ^H? <^> convectively unstable. 

Exponential integrals of the type in (3. 1) have two competing effects. Zeros of the disper- 
sion function ui = ft(k) can lead to exponential growth when 3{f2(fc)} > or cancellation 
and decay due to rapid oscillation when 3{fi(fc)} = for large t. To ascertain whether 
the system is absolutely or convectively unstable one needs to discover which of these op- 
posite tendencies dominates. A number of methods for distinguishing between convective 
and absolute instabilities have been suggested, dating back to the work of Sturrock [12] and- 
Briggs [13]. See also [34, 35, 36]. For completeness we outline the general criteria below. 

Here we assume that T>(uj, k) is known explicitly. The w-integral in (3.1) is along a 
contour that lies above all the zeros of T>(ui, k) for each fixed, real k and we further assume 
that T>(ui, k) is entire in (uj, k) above this contour. Hence, for t > the w-integral may be 
carried out by closing the contour in the lower half-plane and summing over the residues 
of the dispersion function expanded at each of its roots. Assuming the roots of T>(uj, k) are 
simple (multiple roots do note pose a serious difficulty [32]), the resulting integral can be 
written as 

N f°° e i( kx -n n (k)t) 

= (3 - 2) 

where the sum is over the N zeros of the dispersion function 

v(n n (k),k) = o, = . 

The problem is to determine the long time behavior of (3.2) for which the method of 
steepest descent is applicable (cf. [37]). For this, we restrict ourselves to the point moving 
with speed V, x = Vt. Then, by suitable deformation of the real line to the steepest descent 
contour, the dominant contributions arise from the saddle points of the exponent satisfying 

■^r^n(k n ,m) = V , n = l,...N, m = l,...,M n , 

allowing for multiple saddle points along each branch of the dispersion relation. Note that 
the zeros of V, double roots of the dispersion function, do not contribute appreciably to the 
integral because they cancel in the sum (3.2). Using the method of steepest descent, one 
recovers the dominant long time behavior 



G(Vt,t) - 0(e 7 """ t/Vi) , t 
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where 

7max(V) = max3{fi„(fc„ !rn ) - Vk n , m } . 

Thus, if 7 m ax > 0, then an impulse perturbation at t = 0, x = grows without bound along 
the line x = Vt and the instability is absolute. Otherwise, if j max < 0, the perturbation 
decays along the line x = Vt and so the instability is convective. These have been referred 
to as the Bers-Briggs criteria [13, 32]. 

3.2. Simplified criteria for the separatrix of soliton instabilities. Many previous stud- 
ies applied the general criteria for classifying instabilities to dissipative systems (plasma 
physics, viscous fluids, etc.) where the dispersion relation was known explicitly. Given 
explicit (and sufficiently simple) dispersion relations, the analysis of the stationary points can 
be carried out directly. However, the dispersion relation Sl(fc) is unknown for dark solitons 
of the NLS equation. It can be computed numerically, but this makes the analysis of saddle 
points in the complex- k and / or complex-w planes quite challenging. Fortunately, as derived 
below, there are simplified analytic criteria for the transition point between absolute and con- 
vective instabilities of NLS solitons that rely solely on computations of the dispersion relation 
for real k. 

Using the Laplace transform in eq. (2.5), the linearized evolution of an initial L 2 (M. 2 ) 
perturbation tpo(£, y) to the dark soliton satisfies 

V>(t, V,t) = [ e~ iut OC + iuif 1 ¥>o(£, V) du, 

where the Bromwich contour Cb lies above all eigenvalues of C. In order to investigate the 
unstable transverse dynamics in (y, t), we project onto the eigenfunction fo and perform the 
contour integration over Cb resulting in the following representation of the dynamics 

■ roc j(%-n (fc)t) 

The integral is taken over (0, oo) by use of the invariance k — > —k of the eigenpair 

(n (fc), 

By performing a Galilean shift in the NLS equation (1.1) as 

iP(x,y,t)^Tp'(x,y,t) = e i (- w y- wH MTp(x,y + wt,t), (3.3) 
the dispersion relation for transverse perturbations becomes 

fio(fc) -> to(k) = -wk + n (k) , (3.4) 
where — w is the flow speed parallel to the plane of the dark soliton (2.4). This is equivalent 
to investigating the behavior of the perturbation in eq. (3.3) along the line y — wt. With 
this substitution, we consider eq. (3.3) whose long time asymptotic behavior requires the 
evaluation of 

/■^cutoff „ — iQ(k)t 

I ( t )= n ru\ f o(k)dk, (3.5) 

Jo fio(fc) 

where the dependence on ^ is suppressed. The integral over (fc cu toff , oo) is negligible because 
the dispersion relation is purely real (the stationary phase method yields algebraic decay in t, 
cf. [37]). Introducing the change to a complex variable z — Q,(k), eq. (3.5) becomes 



J (*)= / o ( uv? M z ) dz > (3-6) 
where VL'(z) = —w + £l (z) and the contour is C = {z = f2(fc) | k £ [0,fc cuto ff]}- Two 
distinct possibilities arise. 
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FIG. 3.2. Integration contours C (solid curves) in the complex z plane for eq. (3.6) and the real interval 
[— wkcutoff, 0] (dashed lines). The filled circles correspond to poles of the integrand where Q'(k) = f2g(fc) — w = 
0, z = f2(fc), which in (a) prevent the smooth deformation of C to [— tufe cu t ff , 0] giving rise to an absolute 
instability. Parameter values are v = 0.5, ui cr ~ 0.535. (a) w = 0.5 < ui cr . (b) w = 0.6 > ui C r- See also 
Fig. 3.3. 



1. A zero of fl' gives a residue contribution to Cauchy's theorem when deforming C 
to the real interval [— wk cuto g, 0] as in Fig. 3. 2(a). In this case the integral diverges 
exponentially as t —> oo and the instability is absolute. 

2. The zeros of 57' do not lie between C and the real line as in Fig. 3. 2(b) (they may lie 
on the real axis) so that there is a smooth deformation of the contour C to the real 
interval [— wk cu toS, 0], In this case the integral decays to zero as t —> oo and the 
instability is convective. 

0.7 

0.65 



0.5 
0.45 

0.25 0.3 0.35 0.4 0.45 0.5 

k 

FIG. 3.3. Plot ofQ' (k)for real k > fc cuto ff ~ 0.230 and v = 0.5. The minimum of this curve corresponds 
to the coalescence of the poles in Fig. 3.2 and the critical transverse flow speed w CT at which the instability changes 
from absolute to convective. The dashed lines correspond to the values of w used to compute Fig. 3.2(a) (lower, 
absolute instability) and Fig. 3.2(b) (upper, convective instability). 

As discussed in the previous section, the saddle points fi'(fco) = ^o(^o) — w = give 
the long time asymptotic behavior ip ~ O(e % ™ ko ' t / 'Vt), t —> oo. As the transverse flow 
speed w is varied, the type of instability changes from absolute to convective. The transition 
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from absolute to convective instability occurs at w = w CI when two zeros of il'(k) merge 
on the real line. That is, they form a double zero so that f2"(fc cr ) = f2o(fe cr ) ~ an< ^ 
^o(fccr) is minimum. This behavior is depicted in Figs. 3.2 and 3.3 with 3.2(a) showing two 
complex conjugate zeros for w < w CI while 3.2(b) reveals their splitting into two real zeros 
for w > w CI . These real zeros are depicted in Fig. 3.3 for w > w CI . By an appropriate 
choice of the branch cut, one can show that ilo(k) = ft^k*) so that complex zeros of 
fV(fc) = O,' (k) — w come in conjugate pairs. This proves 

PROPOSITION 3.2. The critical wavenumber k cr and critical transverse velocity w CI for 
the transition between absolute and convective instability are real. They satisfy the simplified 
criteria 

-Q^jikcr; w cr ) = , (3.7a) 

QCl 

^(fccnWcr) = 0. (3.7b) 

These conditions were first proposed in [11]. 
Corollary 3.3. 

=0, (3.8a) 



dk 2 



W CT = ^r( fc cr) ■ (3.8b) 



The proof follows from (3.4) and (3.7). 

When the transverse flow speed is subcritical, w < w CI , the dark soliton is absolutely 
unstable and when w > w CI the dark soliton is convectively unstable. The soliton family is 
parametrized by its amplitude v, thus v i-> Wcx(y) forms a separatrix between absolute and 
convective instabilities. The separatrix w cr (f) can also be interpreted as the speed at which 
an initially localized perturbation spreads in time. Thus a convective instability occurs when 
the background flow speed, carrying the perturbation's center of mass, exceeds the speed at 
which the perturbation spreads out. 

In general, the determination of Wc T {v) via (3.8) requires numerical computation. Even 
so, Eqs. (3.8) are much easier to use than the general criteria because the general criteria 
depend on tto(k) over the complex-fc plane whereas (3.8) only depends on Qo(k) for real k. 

3.3. The separatrix in the shallow amplitude regime. The shallow-amplitude asymp- 
totics of the dispersion relation (2. 13) enable us to explicitly compute Qfe and ftkk, determine 
the critical wavenumber k cr and find the separatrix w cr {v) between absolute and convective 
instabilities. Here it is convenient to use the wavenumber scaling (see Appendix A) 

k = v 2 p . 

The asymptotic dispersion relation (3.4) becomes 

n(p; w) ~ -v 2 wp + ^(2V3p - 3) 1 / 2 + u 5 P 2 (^~P] 
= -v 2 wp + v 3 A (p) + ^Axip), 0<^<Cl. 
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The simplified criteria (3.8) give 

fccr = v 2 p CI = is 2 (po + v 2 pi) + 0(v 4 ) , 

n kk (k cr ) ~ ~A£(po) + v [Ai"(po)pi + Ai'(po)] , 

^fe(fc cr ; w cr ) - -w cr + ^Aq(po) + ^ 3 [Ao(p )Pi + AiOo)] ■ 

Equating like coefficients of v and using (3.8), yields 

PROPOSITION 3.4. The first order asymptotic approximation of the critical velocity and 
wavenumber are 

W CI =v+ ^ +0(i/ 5 ), (3.9a) 



KP 



9 

NLS correction 
1 2 4 



V3 3^ 



OfV 5 ) . (3.9b) 



KP NLS correction 

For comparison, Fig. 3.4 shows the numerical solution of the system (3.8) and the asymp- 
totics in (3.9). The numerical details are presented in Sec. 5.1. 

3.4. Convective/absolute instabilities of spatial dark solitons. The natural reference 
frame for studying the convective or absolute nature of soliton instabilities is the one moving 
with the soliton. In this reference frame, both the soliton density and velocity are independent 
of time. The dark soliton is referred to as a spatial dark soliton. Such structures arise, for 
example, in the context of flow past an impurity [9, 38], flow over extended obstacles, and 
dispersive shock waves [39, 40]. 

The spatial dark soliton in (2. 1) satisfies v — 0, which determines the phase jump 
csin/3 - dcos/3 

cos0 = . (3.10) 

Vp 

This soliton is uniquely determined by four parameters rather than five. We use the back- 
ground density p and background velocity (c, d) as three of these parameters along with either 
the normalized soliton amplitude < v = \ sm<p\ < 1 or the soliton angle < f3 < ir/2, the 
two being related via (3.10) through 

v 2 = l _ (csin/3- dcos/3) 2 

P 

The spatial dark soliton exhibits either an absolute or convective instability depending 
on the Mach number of the background flow (2.2) and either the amplitude v or, equivalently, 
the soliton angle (3. By moving in the reference frame £ = x — Kt of the normalized dark 
soliton (2.4), the background flow has velocity — n normal to the soliton and velocity — w 
parallel to the soliton. The critical Mach number of the background flow and its first order 
asymptotic approximation are 



M a (y) = ^J^+w c ,{ V f = v/l-z/ 2 + Wcr H 2 (3.12a) 
(3 =° 1 + \v* + 0(v 6 ) . (3.12b) 

Transverse perturbations are absolutely unstable for M < M CI and convectively unstable for 

M > M ct . 

We also compute the critical Mach number's dependence on the soliton angle by use of 
eq. (3.11) leading to a transformation between v and f3 

w 

tanp = ; . 
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FIG. 3.4. (a) The separatrix w CT (u) between absolute and convective instabilities. For speeds w > w CT 
(w < w cr ), the dark soliton is correctively (absolutely) unstable, (c) The critical wavenumber k cr satisfying the 
condition (3.7a). (b) and (d) are the differences in the numerically computed values ofw cr , k CI and the first order 
approximations in (3.9a) and (3.9b), respectively, showing the expected error settlings. 
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1.44 
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1.44 



computed 
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first order 
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FIG. 3.5. Critical Mach number M cr (Eq. (3.12a)) and its asymptotic approximations (M CT = 1 for KP, 
Eq. (3.12b) for the first order correction) as functions of the soliton amplitude u (a) and the spatial soliton angle (5 
(b). For M > M cr , the spatial dark soliton is convectively unstable, otherwise it is absolutely unstable. 



Therefore, shallow spatial dark solitons at critical velocity w — w cr (v) have a small j3 
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v + (D(v 3 ) < 1, so that 

M cr (/3) = l + ^/3 4 + 0(/3 6 ). (3.13) 

Figure 3.5 shows the numerically calculated dependence of M cr on v and f3 and compar- 
isons with the asymptotic results (3. 12b) and (3.13). Combining the asymptotic result (3.12b) 
with these computations leads to 

CONCLUSION 3.5. The transition between convective and absolute instability for spatial 
dark solitons always occurs at supersonic speeds M CI > 1. A sufficient condition for a spatial 
dark soliton with background Mach number M to be absolutely unstable is 

M < 1 . 

A sufficient condition for a spatial dark soliton with background Mach number M to be 
convectively unstable is 

M > M c ,{v = 1) « 1.4374 . 
Additionally, v t— > M CT (y) is monotonically increasing. In sum, 

1 < M CT < 1.4374 . (3.14) 

REMARK 3.6. In [11], the bounds 1 < M CT < 1.46 were obtained. The leading order 
termin eq. (2.13)was used to show that M CT ~ 1 in the shallow regime. Equation (3. 12b) im- 
proves the lower bound on M CT and demonstrates that M cr is strictly supersonic for all finite 
soliton amplitudes. The upper bound 1.46 in [11] was calculated from a rational approxi- 
mation of the spectrum for large soliton amplitudes [29]. Equation (3.14) gives the accurate 
upper bound. 

4. Oblique dispersive shock waves. In a dispersive fluid where dissipation is negligi- 
ble, a jump in the density /velocity may be resolved by an expanding oscillatory region called 
a dispersive shock wave. The Whitham averaging technique [41] has been successfully used 
to describe a DSW's long time asymptotic behavior in a number of physical systems, for 
example [42, 43, 17, 44, 45]. We briefly recap the rudiments of DSW theory. A DSW is 
a modulated wavetrain composed of a large amplitude, soliton edge and a small amplitude, 
oscillatory edge, each moving with different speeds. In the relatively simple case where 
a DSW connects two constant states, the speeds associated with each edge are determined 
by jump conditions [46], in analogy with the Rankine-Hugoniot jump conditions of classi- 
cal, viscous gas dynamics. The jump conditions result from a simple wave solution of the 
Whitham modulation equations connecting the zero wavenumber, soliton edge to the zero 
amplitude, oscillatory edge. The existence of a DSW for a particular jump in the fluid vari- 
ables is guaranteed when an appropriate entropy condition is satisfied. For a left-going DSW, 
we define the leading (trailing) edge to be the leftmost (rightmost) edge - and vice versa for a 
right-going DSW. The sign of the dispersion determines the locations of the soliton and small 
amplitude edges. For systems with positive dispersion such as the NLS eq. (1.1), the soliton 
is a depression wave that resides at the trailing edge of the DSW. 

While DSWs in (l+l)-dimensions have been well-studied, the theory of supersonic dis- 
persive fluid dynamics in multiple spatial dimensions is in its infancy. Perhaps the simplest 
DSW in multiple dimensions is an oblique DSW, which has been studied in the stationary 
[47, 48, 49, 39] and non-stationary [40] regimes (see Figs. 4.1 and 4.4). In this section, 
the analysis from the previous section is applied to the stationary and nonstationary oblique 
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DSW soliton trailing edge to determine the separatrix between convective and absolute in- 
stabilities. In addition, in the weak shock and hypersonic regimes, we find that the jump 
conditions for stationary and nonstationary oblique DSWs are the same. As in classical gas 
dynamics, oblique DSWs can serve as building blocks for more complicated boundary value 
problems. Therefore, understanding the instability properties of oblique DSWs is important 
and relevant to supersonic dispersive flows. This has been further demonstrated by recent 
numerical simulations of NLS supersonic flow past a corner [39, 40]. 

In Sec. 4.1, the jump conditions and instability properties of nonstationary oblique DSWs 
are presented. The following Sec. 4.2 contains a derivation of a stationary oblique DSW in 
the shallow regime, its stability, and comparisons with numerical simulation. Finally, Sec. 4.3 
demonstrates the connections between stationary and nonstationary oblique DSWs. 



v 




x 



FIG. 4.1. Schematic of an oblique DSW. 



4.1. Nonstationary Oblique DSWs. In this section, we first recap the derivation of a 
nonstationary oblique DSW [40] and then discuss its instability properties. 

A schematic of a non-stationary oblique DSW at a specific time in its evolution is de- 
picted in Fig. 4.1. An incoming upstream, supersonic flow is turned through the oblique DSW 
by the deflection angle 9. To accommodate the deflection, the oblique DSW expands along 
the wave angle j3. The leading edge consists of small amplitude waves propagating into the 
upstream flow while the trailing edge is composed of a dark soliton whose amplitude and 
speed are asymptotically calculated from the oblique DSW jump conditions. 

The nonstationary oblique DSW results from the long time evolution of an initial jump 
in the density and velocity component normal to the DSW wave angle /3, in the direction 
hp = (sin/3, — cos/3), and continuity of the velocity parallel to f3, in the direction pp — 
(cos (3, sin /3). We consider the upstream state 

lim p = pi , lim u = (u\, 0) , 

x— > — oo x~ > — oo 

and the downstream state 

lim p — p2 , lim u — (v,2 cos 9, u-i sin 9) . 

x— >+oo x— >+oo 

The normal 1-DSW associated with the dispersionless characteristic Ai = u — ^fp (left-going 
wave) satisfies the simple wave condition [43] 

hp ■ (ui — u 2 cosd, ~u 2 sin#) = 2( v / pJ — ^fpl) ■ (4.1) 
A NLS governed fluid experiences potential flow (see eq. (1.2)). By restricting the spatial 
variation of the solution to the direction hp and integrating the irrotationality constraint v x — 
u y along the direction pp, we obtain the continuity of the parallel velocity component 

pp ■ (ui — U2 cos 9, —u 2 sin 9) = . (4.2) 
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Choosing the reference frame in which the soliton trailing edge is fixed, the speed of the 
soliton edge satisfies [43] 

VK,0)- A /^ = 0. (4.3) 

V pi 

The jump conditions (4.1), (4.2), and (4.3) for the oblique DSW relate the upstream quantities 
pi, Ui and one of the angles 9 or /3 to the downstream quantities p 2 , u 2 and the other angle. 
Introducing the Mach numbers Mj = Uj/^/pj, j — 1,2 along with some manipulation, the 
jump conditions become [40] 

2 

tan(/3 - 6) = — sec /3 - tan /3 , (4.4a) 
Mi 



M2 . . ^Ejggg , (4 . 4b) 

cos(p — y) Mi sin /? 

/92 = PiMi 2 sin 2 /3. (4.4c) 

Further manipulations lead to the equivalent relations 

sin(2/3 - 6>) 2 M x cos(2/3) + 2 sin /3 

cos^ — 



cos(/3 - 6») Mi ' ' y/A + Mf -4Misin/3 ' 

The associated entropy condition is p 2 > Pi, which, when incorporated into the jump 
conditions, gives 

Mi > 1 , < 9 < it , sin" 1 -±-<0<l. 

These state that the upstream flow must be supersonic, the flow always turns into the DSW, 
and the wave angle is larger than the Mach angle sin _1 (l/Afi). The Mach angle is half 
the opening angle of the Mach cone inside of which infinitesimally small disturbances are 
confined to propagate in dispersionless supersonic flow. A convenient way to visualize these 
results is by the M-9-/3 diagram in Fig. 4.2 that relates the deflection and wave angles for 
a given upstream Mach number Mi. Figure 4.2 includes the sonic curve M 2 = 1 (to the 
right/left the flow is sub/supersonic). 

A natural question is whether oblique DSWs with supersonic downstream flow conditions 
are convectively or absolutely unstable. To address this question, we use: 

DEFINITION 4. 1 . Transverse perturbations to the nonstationary and stationary oblique 
DSW are convectively (absolutely) unstable whenever the trailing, dark soliton edge is con- 
vectively (absolutely) unstable. 

See further discussion in Sec. 6. 

Spatial dark solitons exhibit the constraint (3.10). When applied to the oblique DSW 
trailing edge in Fig. 4.1 with background flow parameters (c,d) = y/pM2(cos6,sm6), we 
find 

cos 4> = M 2 sin(^ - 6) . 

Using the jump conditions in eqs. (4.4), we determine the normalized soliton amplitude 

i/(Mi ,d) = sm<p= — . , (4.5) 

Mi sin p 

where f3 is related to 6 by (4.4a). The Mach number of the downstream flow adjacent to the 
soliton is M 2 so the absolute/convective stability criterion (3. 12a) determines the separatrix 

M 2 (Mi, 0) = M cr (» = Vl-^+^M 2 , (4.6) 
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wave angle j3 [deg] 



FIG. 4.2. The M-8-/3 diagram for non- stationary oblique DSWs of the NLS equation (I.I). Each upstream 
Mach number M\ leads to a relationship between 9 and f3. The separatrix curve ( solid) between convectively and 
absolutely unstable solitons is supersonic, i.e. in the M2 > 1 region (left of the sonic line, dashes). The separatrix 
curve asymptotes to the sonic line as (3 — > 90°. 



with v given in (4.5). Conclusion 3.5 implies 

COROLLARY 4.2. Nonstationary oblique DSWs with subsonic downstream flow are 
absolutely unstable. Supersonic downstream flow can be either convectively or absolutely 
unstable. 

This conclusion can also be gleaned from Fig. 4.2. To the right of the separatrix, the 
trailing edge oblique soliton is absolutely unstable because M2 < M cr while to its left, the 
soliton is convectively unstable. The region to the right of the separatrix and to the left of 
the sonic line represents absolutely unstable oblique DSWs with supersonic downstream flow 
conditions. Below we derive additional properties of the separatrix. 

From Fig. 4.2, we observe a minimum wave angle j3 CI , below which the oblique DSW is 
convectively unstable. Setting Mi = M cr in eq. (4.4b) and solving for j3 we find 



_ -2+V4+(4 + Mf)A/ c 2 r 



which has a minimum for Mi = 2y/l + M*, M cr = M ct (v = 1) w 1.4374. We therefore 
have a sufficient condition for the oblique DSW trailing edge to be convectively unstable 



/3</?c 



(i + M cr (i) 2 y 1/2 



34.83°. 



The nonstationary oblique DSW is uniquely determined by the parameters Mi, 6, and 
p\. Thus, given Mi > 1 and < 9 < n, the absolute or convective instability of the 
corresponding oblique DSW's trailing edge is determined by the location of (Mi, 9) relative 
to the separatrix condition (4.6) in the Mi-6* plane as shown in Fig. 4.3. The parameter p\ 
does not affect the absolute or convective nature of the instability. 
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[deg] 

FIG. 4.3. Separatrix curve (solid) dividing the M\-9 plane for nonstationary oblique DSWs into two regions: 
the absolute/convective instability of the trailing edge dark solium. The sonic curve M2(Mi, 9) = 1 (dashed) is 
included for comparison. 



Using the small amplitude result (3. 12a), M CT ~ 1 + |^ 4 , assuming near sonic upstream 
flow Mi — 1 + e, < e = 0(v 2 ) -C 1, and expanding (3, 9, Mi, and M cr , we compute the 
critical angle 

M 2 (l + e, 9 CT ) = M CT + 0(e 3 ) = M CI + 0{v 6 ) , 

4 Wl-^, 0<e«l. 



3^ V 27 

For 6 < 6 CI , the trailing edge dark soliton is convectively unstable and absolutely unstable 
otherwise. Similarly, the sonic angle satisfies 

M 2 (l + e, 9 sonic ) = 1 + 0(e 3 ) = 1 + 0{is 6 ) , 

4 =e 3/2 f 1 - ) , 0<e<l. 



v some ^ ^ 3 

For < S onic> me downstream flow is supersonic and subsonic when 9 > SO mc- F° r the 
narrow window of deflection angles 9 CI < 9 < 9 son [ C , the flow is supersonic and absolutely 
unstable. 

4.2. Spatial Oblique DSWs. We have so far focused on nonstationary oblique DSWs. 
In this section, we construct stationary or spatial oblique DSWs in the weakly nonlinear 
regime (see Fig. 4.4), study their instability properties, and perform numerical simulations. 
This discussion for the NLS equation (1.1) with positive dispersion parallels the developments 
in [47, 48] applied to ion-acoustic waves in plasma, a system with negative dispersion. 

Equations (1.3a), (1.3b), and the irrotationality constraint due to potential flow are con- 
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sidered in the (2+0)-dimensional case 

{pu) x + (pv) y = , (4.7a) 

uu x+ vu y + p x = \(^^-^-fi^j , (4.7b) 
v x — u y = . (4.7c) 

We seek a special class of solutions that are related to supersonic flow past a sharp corner 
or wedge. For this, we treat y as a time-like variable and consider the "initial conditions" at 

y = o 

p(x,0) = ( X< ° n , (4.8a) 

I n) f Ml, X < ,. „, . 

u(x,0) = i ,— n n ? (4.8b) 

1 y/p!M 2 cos 0, £ > ' 



w(a;,0) = <j _°' . . 2;< !1 • (4.8c) 
> VpiM 2 sin 0, x>0 v y 



The well-posedness of this initial value problem is plausible in the supersonic regime, Mj > 
1, j = 1,2, due to the hyperbolicity of the dispersionless equations (see e.g. [27]). We 
seek a stationary, oblique DSW solution in the supersonic and weakly nonlinear regime < 
pi — 1 -C 1. For this, we apply the method of multiple scales 

p = 1 + ep (1) + e 2 p {2) + ■■■ , (4.9a) 
u = Mi - eu {1) + eV 2) + • • • , (4.9b) 
v = sv {1) + e 2 v {2) + • • • , (4.9c) 

in the transformed variables 

f = e 1/2 [x~ (A/ 2 - l) 1/2 y], T = e 3/2 y. (4.10) 
This particular choice is motivated by the line £ = const whose angle with the x axis is 
the Mach angle sin _1 (l/Ajfi) for small amplitude wave propagation in the upstream flow. 
Equating like powers of e leads to 



f ) +M 1 p^-{M 2 -l)i V ^ = 0, 
0(si): _ Mlu m +/9 W = o, 

«J X) - (M? - l)3uf } = 0. 

The solution incorporating the initial conditions (4.8) is 

p« = M lU « w« = (Af? - l)*uW , (4.11) 
with uW determined at the next order: 



0(e*) 



4 2) + M lP f ) - (M 2 - 1)^ 2) - (p(%«) ? 

+4 1) -(M 2 -i)3(p( 1 V 1 )) e = o, 

MlU f + pf + uWu^ + (M 2 - l)^(i)4 x) = ip« 

wf ) + (M 1 2 -l)34 2) +4 1) = 0. 
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Inserting eqs. (4.1 1) we obtain the KdV equation 

U (D ^L_ U (D U W + = o . (4.12) 

2(M 1 2 -l)i ? 8(Mf-l)i f " 
It is convenient to consider the transformed variables U, C as 

U = - 3M * rU W+l, Z= M fo- T \. (4.13) 
(Af?-1)* 2(M 1 a -l)* 
Then, eq. (4.12) becomes the KdV equation with negative dispersion 

?7 T + C/C/ c + U ac = . 

The initial data in (4.8) maps to the Riemann problem 



u(C,o) = 



i c < o 
o c > o 



This dispersive Riemann problem was solved by Gurevich and Pitaevskii in 1974 [42]. The 
result is a DSW with the trailing edge, small amplitude wave speed ct = — 1 and leading 
edge, soliton speed cl = 2/3. The leading edge soli ton amplitude is 2 corresponding to the 
KdV soliton speed/amplitude relation. The oscillatory part of the DSW, for r sufficiently 
large, has the approximate form [42, 15] 



a ,'K[m(C/r)} 



U(C, t) ~ m(C/r) - 1 + 2dn z ^ 1 ^ " <KC, r);m( C /r)j , r » 1 , (4.14) 

where dn is a Jacobi elliptic function and A" [to] is the complete elliptic integral of the first 
kind. The elliptic parameter to(£/t) is the self-similar, simple wave solution to the Whitham 
modulation equations given implicitly by 

C I, . 2 (l-m)K\m\ 
- = -(1 + to)--to- v 



r 3 V y 3 E[m] - (1 - m)AT[m] 

where B[m] is the complete elliptic integral of the second kind. The phase is determined 
through 



2 / 3 dz 



C/r K[m{z)] ' 

To obtain the NLS oblique DSW solution in its unsealed form, we use the transformations 
(4.13), (4.10) along with the substitutions (4. 1 1) to match the asymptotic solution (4.9) to the 
initial conditions (4.8). The deflection angle is related to the small parameter e via 

(M?-l)l 

1 J «1, (4.15) 

so that weak spatial DSWs correspond to a small DSW deflection angle. Then the relationship 
between the downstream and upstream variables takes the asymptotic form 

P2-1+ , 0, (4.16a) 

(Ml - 1)5 

/ m 2 I o \ 
M a ~Afi 1 ^ — -6\, 0<6»<1. (4.16b) 

The KdV DSW speeds ct = — 1 and cl = 2/3 correspond to the slopes of the oscillatory 
region's boundaries which we transform to the leading and trailing angles f3 + , /?_, respec- 
tively, for the stationary oblique DSW. Using the transformations (4.10), (4.13), and (4.15), 
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the oblique DSW angles take the asymptotic forms 



/3_ ~ sin" 1 — 



P+ ~ sin" 1 — 





Ml 


Mi) 


+ 2(M 1 2 - 


±) 


3M? 


Mi J 


+ Ml-l 



e , o < <?« 1 



(4.17a) 
(4.17b) 



Finally, the trailing edge soliton amplitude and phase jump 2<fr with the angle (3 have the 
asymptotic form 



p 2 - p(x, xt&nP ) = v /7jjsin< 



(4.18) 



This DSW solution is plotted in Fig. 4.4 and approximates a stationary, weak, oblique DSW 
forNLS. 




FIG. 4.4. Example spatial oblique DSW with small deflection angle 8 constructed from the asymptotic solution 
(4.14). 

Equations (4.16) and (4.17) are the jump conditions for weak, stationary oblique DSWs. 
These can be used to approximately solve the problem of supersonic flow over a corner with 
angle < 8 <C 1. Additionally, due to symmetry arguments, two stationary oblique DSWs 
approximately solve supersonic flow over a wedge as in [39]. Figure 4.5 shows the numerical 
solution of eq. ( 1 . 1) for supersonic M\ = 2 flow past a corner with angle 9 = 9° after the flow 
pattern has reached a quasi-steady state (see Sec. 5.2 for the numerical details). Sufficiently 
close to the corner, the structure of the numerical solution resembles the asymptotic oblique 
DSW shown in Fig. 4.4. Further away from the corner, the first sign of instability occurs 
along the trailing, dark soliton edge leading to the generation of vortices. This provides 
some justification for our definition of oblique DSW instability in Def. 4.1. Furthermore, 
we observe that the vortices are convected further away from the corner as time progresses 1 . 
In a previous work [40], the authors performed numerical simulations of NLS supersonic 
flow past a corner for a large number of flow configurations, observing similar, stable pattern 
formation in some cases. Flow configurations where the instability overwhelms any stable 



The vortex pattern eventually stabilizes at a fixed distance from the corner. A recent study [50] of NLS dark 
soliton convective/absolute instabilities has some independent results that overlap with ours in Sections 3.2 and 3.4. 
This work also gives a further description of perturbation convection along the soliton. The effective group velocity 
of the perturbation along the soliton is found to be equal to the critical flow speed (here ^/p^Wc?). However, 
convective instability theory does not explain the numerically observed stabilization of vortex formation at a fixed 
distance from the comer. 
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M x = 2 


9 


£ 


M 2 


Pi 


13- 


P + 


theory 


3° 


0.6 


1.82 


1.12 


32° 


42° 


numerics 


3° 


0.6 


1.84 


1.12 


32° 


39° 


theory 


6° 


1.3 


1.64 


1.24 


34° 


54° 


numerics 


6° 


1.3 


1.67 


1.26 


34° 


49° 


theory 


9° 


1.9 


1.46 


1.36 


36° 


66° 


numerics 


9° 


1.9 


1.51 


1.40 


37° 


62° 



Table 4.1 

Comparison between the asymptotic results ofeqs. (4.16) and (4.17) and numerical simulation of supersonic 
flow over a corner. 



pattern formation were also observed. We identify these two flow regimes with convective 
and absolute instability of the oblique DSW. 




X 



FIG. 4.5. Numerical simulation of supersonic flow past a corner with 9 = 9°, Mi = 2, at t = 400. The 
color scale is chosen to visually resolve the small amplitude oscillations. 

Table 4.1 summarizes the asymptotic estimates in eqs. (4.16) and (4.17) compared with 
the numerical computations showing excellent agreement, even for fairly large corner angles 
and when the "small" parameter e is larger than one. 

The trailing edge dark soliton is shallow. Therefore, using the theory developed in 
Sec. 3.3 and eqs. (4.16b), (3.12b), the oblique DSW is convectively unstable when 

M 2 > Af cr (» or M x [l + 0(61)] > 1 + 0(9 4 ) , < 9 < 1 , 

because v = sin</> ~ 0(9) from eq. (4.18). As long as Mi > 1, independent of the corner 
angle 8, using Conclusion 3.5 gives 

COROLLARY 4.3. For NLS supersonic upstream flow M\ > 1 and sufficiently small 
corner angles < 9 <C 1, the oblique DSW emanating from a sharp corner is convectively 
unstable. 



24 



M. A. HOEFER AND B. ILAN 



4.3. Relationship between stationary and nonstationary oblique DSWs. As shown 
in the previous section, stationary oblique DSWs can be physically realized as the solution 
of a two-dimensional boundary value problem involving supersonic flow. In contrast, the 
nonstationary oblique DSW studied in Sec. 4. 1 results from the solution of an initial value 
problem. As we now demonstrate, the downstream flow conditions for the stationary and 
nonstationary oblique DSW are the same in two asymptotic regimes: weak shocks and hy- 
personic flow. 

The downstream flow conditions and the stationary trailing edge soliton in both the sta- 
tionary and nonstationary oblique DSWs are characterized by the deflection angle 6, the 
wave angle /3~ or f3 for the nonstationary case, the Mach number A/2, and the density p^. 
These properties are related via the oblique DSW jump conditions. For weak oblique DSWs, 
we assume a fixed upstream supersonic Mach number M\ > 1 and small deflection angle 
< 9 <C 1 as in Sec. 4.2. By a standard asymptotic calculation, an expansion of the jump 
conditions for the nonstationary oblique DSW in eqs. (4.4) in the form 

pi = 1 , P 2 = 1 + 0(9) , M 2 = M x + 0{9) , (3 = sin- 1 l/M x + 0{9) , 

gives precisely the same result as that obtained for the stationary oblique DSW in eqs. (4. 16a), 
(4.16b), and (4. 17a). 

The hypersonic regime assumes the large Mach number scaling Mi ^ 1 and small 
deflection angle 9 — 0[\jM\). In this asymptotic regime, the jump conditions (4.4) become 



P»=(j^ L + 1 ) (4.19a) 
2Mi 

M2 =2Ti + °W> (4 - 19b) 

/3 = ~ + -3- + Oil/Ml) , Mi>l, 0<9 = O(l/M 1 ), (4.19c) 
i Mi 

where we have assumed that pi = 1. In [49, 39], stationary oblique DSW solutions of 
eqs. (4.7) and (4.8) in the hypersonic regime were constructed asymptotically. The classical 
notion of hypersonic similitude [51] applies so that the (2+0)-dimensional stationary problem 
was asymptotically mapped to a (l+l)-dimensional problem for the NLS equation. Station- 
ary, supersonic flow past an extended obstacle is then related to a piston problem, which can 
be solved analytically in the case of a sharp comer (constant piston speed) [52] and for more 
general profiles [39, 53]. The results for the stationary oblique DSW are the same as those 
computed asymptotically for the nonstationary case in eqs. (4.19) when M X Q < 2. The case 
M\8 > 2 corresponds to a novel feature of the dispersive piston problem where the oblique 
DSW experiences cavitation and the DSW forms an oscillatory wake [52], not captured by 
the jump conditions (4.4). Combining these results with Corollaries 4.2 and 4.3 demonstrates 
CONCLUSION 4.4. For weak (small deflection angle < 6 <C 1, fixed upstream Mach 
number Mi) or hypersonic (Mi > 1, 9 = 0(1/Mi), M x 6 < 2) oblique DSWs, the non- 
stationary and stationary flows have the same asymptotic downstream flow properties and 
trailing edge soliton amplitudes/angles. In these regimes, the oblique DSWs are convectively 
unstable. 

5. Computational techniques. In this section, we present details of our numerical 
methods for computing the spectrum of transversely unstable perturbations as well as deriva- 
tives of the dispersion relation via adjoint methods (Sec. 5.1). Direct numerical simulations 
of NLS supersonic flow over a corner are explained in Sec. 5.2. 
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5.1. Computing the spectrum and derivatives of the dispersion relation. Accurately 
computing the spectrum of the linearized NLS equation (1.1), finding the maximal growth 
wavenumber (2.14) and the critical wavenumber (3.7a) require a fine grid and sufficiently 
large computational domain. This turns out to be challenging, especially in the small am- 
plitude regime. To achieve this, we employ a combination of computational and analytical 
techniques explained below. 

• The linearized operator in (2.7) is realized using the centered, fourth order (sparse) 
finite difference stencil in £ for the Laplacian and other derivative operators. Zero 
Dirichlet boundary conditions are embedded into the associated matrix. We find that 
a domain size of 11 /v serves well (increasing the domain size has negligible effect 
on the results). 

• For accuracy, the number of grid points along the transverse direction £ should scale 
as Ijv. As v decreases from 1 to 0.01, we use 2 9 — 2 13 grid points. Using fewer 
points can lead to completely wrong results, either because k niax — > or because 

fccr -> fccutoff+ as V ->• 0. 

• The discrete eigenvalue fio (k) and its associated localized eigenfunction fo (£; k) are 
computed using Matlab's sparse eigenvalue solver ( ' eigs ' with ' SM' ). 

• One approach is to compute fto(k) on a grid of k values (as for Fig. 2.1). Then, 
Q'o(k) (resp. ft'^k)) can be computed using finite differences and minimized on 
the k grid to find k max (resp. fc cr ). This method turns out to be computationally 
expensive. To overcome these challenges, an accurate and fast method is explained 
below. 

Recall the eigenvalue problem (2.7). As discussed previously the discrete spectrum of 
JL consists of two simple eigenvalues of opposite signs, ±Slo(fc) (we choose the positive 
sign), with the associated eigenfunction fo(£; k). Our main goal is to compute fc max such 
that ft' Q (k max ) = 0, and k cr such that (fccr ) = 0. This is achieved using the following 
algorithm: 

1. Compute the discrete spectrum at some (initial) k, i.e. tto(k) and fo(£; k). 

2. Apply adjoint methods to find exact expressions for Og(fc) and fi '(fc), i.e. Eqs. 
(5.5)-(5.6) below. 

3. Repeat steps 1-2 using a root finder to converge to fc max and k CI . 

4. Compute r max = 3{r2o(fcma X )} and/or£l CI = il (k CI ) and w CI = n' Q (k cr ). 

We proceed to derive the relevant expressions. Use will be made of the standard Pauli 
matrices, 





"0 


1 




'0 -{ 




'1 " 




1 





) CT 2 = 


i 


, ^3 


-1 



(5.1) 



the reflection operator R, Rg(x) = g(—x), and the adjoint of an operator will be denoted by 
(•) f - 

Differentiating (2.7) with respect to k gives 

(a 2 L + h 2 <j 2 + n ^) f£ = -(ka 2 + (l' )f , (5.2) 

where (■) denotes differentiation with respect to k. Solvability requires that (ka 2 + f^o)fo be 
orthogonal to the nullspace of the adjoint operator to the left-hand side of (5.2). In Appendix 
C we prove that this nullspace can be characterized as follows. 

LEMMA 5.1. For k g C \ {0, ±fc cuto ff } and (flo(k), fo(£; k)) an eigenpair satisfying 
(a 2 L + \k 2 a 2 + fi ^ f = , (5.3) 
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we have 



1 



ker | [a 2 L Q + ^k 2 a 2 + n Q j j = span{i?cjif *} . (5.4) 

Using Lemma 5.1 and taking the inner product of (5.2) with fieri fg, the solvability con- 
dition reads 

Differentiating (5.2) with respect to k gives 

<J2L + \k 2 a 2 + fi ) fo - -(fc<7 2 + %K - (*s + fio)fc • 

Using the solvability condition we conclude that 

o( j ~ (fo.^xfo*) (fo,^^) ■ (5 - 6) 

In summary, we compute fid(fe) and Oq (/s) using Eqs. (5.2), (5.3), (5.5), and (5.6). These 
computations are accurate and fast. The most time consuming operation is the computation 
of the discrete spectrum of L . 

5.2. Numerical solution of the NLS equation. In Section 4.2 we presented the nu- 
merical solution of supersonic NLS flow past a corner. The technique used was the same 
as that presented in [40]. We introduce a linear potential with large contrast that acts as a 
penalization to flow outside the domain. Such volume penalization methods are well-known 
in classical fluid dynamics (see, e.g., [54]). In the context of BEC and optics, superfluid 
flow around obstacles or boundaries are realized in practice using electromagnetic waves or 
a variable refractive index, both modeled as a spatially varying, linear potential. The benefits 
of this numerical technique include the use of a regular, Cartesian mesh and highly accurate 
pseudospectral derivative calculations. 

The time-dependent NLS / GP equation (1.1) with a linear potential 

ii/H = ~\{i>xx + i> yy ) + V(x, y, t)i> + iVfV* , (5-7) 
was solved numerically using a pseudospectral, Fourier spatial discretization and a fourth or- 
der Runge-Kutta explicit time stepper. These computations were performed on a rectangular 
mesh of N x N y equispaced grid points within the domain [—L x , L x ] x \—L y , L y ]. Our choice 
of the potential 



V(x,y,t) 



1 - H^(L X - |z| - 6)H^(L y - \y\ - 6) (5.8a) 



H^y - C(x - Mit)) 



(5.8b) 



C(0 = - tan(0) [H^xt - £) - H„(l)\ - L y + 8 , (5.8c) 

^(0 = ^ + ^tanh(£//i), (5.8d) 

models the boundary conditions corresponding to flow over a corner and also serves to "lo- 
calize" the solution so that a pseudospectral, Fourier discretization with periodic boundary 
conditions can be employed. An example potential is shown in Fig. 5.1. The time-dependent 
potential corresponds to a moving ramp. The function is a regularized Heaviside step 
function with transition width /1. The terms on line (5.8a) effect the localization of tp to 
within S of the domain boundaries. The terms on lines (5.8b) and (5.8c) correspond to a 



DARK SOLITONS, DSWS, AND TRANSVERSE INSTABILITIES 



27 



moving ramp with corner angle 8 and apex located at (xq — Mjt, —L y + S). When the sec- 
ond corner at x — x — Mjt + I is reached, the ramp flattens and continues as a straight 
line. The initial condition is the nonlinear, stationary ground state of eq. (5.7) with poten- 
tial V(x, y, 0) computed by the spectral renormalization technique [55] with the unit density 
constraint \ip(0, 0, 0)| 2 = 1. The potential contrast Vq is taken sufficiently large so that the 
density is effectively zero where V(x, y, t) « Vq. Time integration was carried out until the 
corner reached the left boundary. Near the corner, the flow approximates a "pure" oblique 
DSW as shown in Fig. 4.5. Parameter values for table 4.1 are N x = 4000, N y = 1000, 
L x = 2000, L y = 500, V = 20, 5 = 2, x = 1000, I = 1000, fi = 2, and a time step of 
0.05. The simulation depicted in Fig. 4.5 results from N x = 3200, N y = 1600, L x = 800, 
L y = 400, V = 20, 6 = 2, x a = 400, 1 = 400, n = 2, and a time step of 0.01. 




FIG. 5.1. An example potential V(x,y,t) with Vo S> 1 used to model numerical simulation of supersonic 
flow past a corner. Regions where V(x, y,t) is large correspond to negligible density. The ramp moves to the left 
with speed Mi leading to oblique DSW formation. 

6. Discussion and conclusions. One of the motivating questions for this study was the 
nature of convective versus absolute instabilities of dark solitons. In general, the character- 
ization of the instability type requires knowledge of the dispersion relation for a range of 
wavenumbers in the complex plane. Unfortunately, the exact discrete spectrum (and hence 
dispersion relation) for NLS dark solitons is unknown. The formal analysis presented in [1 1] 
led to greatly simplified criteria for determining the instability type, which involve only the 
imaginary (stable) portion of the spectrum. 

In this study, the underlying assumptions behind the simplified criteria are exposed and 
justified using a combination of rigorous results (Theorem 2.1 and Lemma 5.1), shallow 
amplitude asymptotics (Prop. 2.3), and computations of the spectrum. Consequences of the 
small-amplitude asymptotics and numerical computations are the first order corrections to 
the maximal growth rate and associated wavenumber (Corollary 2.4) and dependence of the 
critical Mach number on the soliton amplitude (Conclusion 3.5). Applying Conclusion 3.5 
to the soliton trailing edge of oblique DSWs, we conclude that subsonic oblique DSWs 
are always absolutely unstable, whereas supersonic oblique DSWs can be absolutely 
or convectively unstable (Corollaries 4.2 and 4.3). In addition, the relationship between 
stationary DSWs (corner B VPs) and nonstationary DSWs (Riemann IVPs) is studied. In both 
cases, the DSWs are found to have the same downstream flow properties in the shallow and 
hypersonic regimes (Conclusion 4.4). 

It is worth contrasting these results with oblique shock waves in classical gas dynam- 
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ics. Supersonic classical shock fronts in gas dynamics are linearly stable when they satisfy 
the Lax entropy condition [23]. For the boundary value problem of supersonic flow past a 
sharp corner, the oblique shock is stable if and only if the downstream flow is supersonic 
[24, 25]. As far as we know, the distinction between absolute and convective instabilities 
in the subsonic case has not been elucidated. We note that recent experiments in another 
viscous medium (granular material) exhibit the stable excitation of oblique DSWs with both 
supersonic and subsonic downstream flow conditions [56]. 

Several questions and open problems related to this study are mentioned below. 

The nonstationary oblique DS W consists of a slowly modulated elliptic function solution 
to NLS. How to study the stability or instability of this coherent structure is not immediately 
obvious given its expanding nature and asymptotic representation. The notion of instability 
we consider here is centered upon the properties of perturbations to the stationary, trailing 
dark soliton edge. This is a natural criterion because the soliton trailing edge corresponds to 
the largest oscillation in the DSW, hence nonlinear effects are strongest there. Another moti- 
vation for this choice comes from the numerical simulation of supersonic flow over a corner 
where the instability first appears along the trailing edge soliton. However, to gain a more 
complete understanding of DSW instabilities, one should develop an analysis of absolute and 
convective transverse instabilities of elliptic function solutions. This suggests the more gen- 
eral study of convective/absolute instability for systems with continuous bands of unstable 
modes. We are not aware of any previous work in this direction. 

Careful computations of the spectrum suggest that Conjecture 2.2 is true. However, to 
the best of our knowledge, it has not been proven rigorously. It may be possible to do this by 
reducing the problem to an ODE, where Sturm-Liouville theory is applicable (cf. [57, 30]). 

It would be interesting to extend these results to systems with negative dispersion such 
as shallow water waves where the KP-II equation is valid in the small amplitude regime. In 
contrast to KP-I studied here, line soli tons are linearly stable [28]. Are oblique DSWs in 
systems with negative dispersion stable? 

Acknowledgments. The authors thank Mark Ablowitz for inspiring remarks. The au- 
thors also thank Anatoly Kamchatnov for constructive discussions and sharing his recent 
manuscript [50]. 

Appendix A. Eigenvalue asymptotics. 

We seek the dispersion relation Q(k; v) of Eq. (2.7) for unstable transverse perturbations 
to the shallow (0 < v -C 1) dark line soliton (2.4). Rather than perform asymptotics directly 
on (2.7) it is convenient to consider the eigenvalue problem in fluid variables (1.2). The 
soliton solution (2.4) takes the form 

p(x, y, t) = p.(0 = 1 - ^ 2 sech 2 (C) , (A.la) 

u(x,y,t) = u s (C) = o o o i (A.lb) 

^ ' ^' sinh 2 (C) + ^cosh 2 (C) 

v(x, y,t) = , £ = v{x — Kt) . (A.lc) 

Applying multiple scales to ( 1 .3) leads to the KP-I equation for weakly nonlinear excitations 
of (1.1) to the uniform state p = 1 (cf. [58]). The scalings involved motivate the following 
representation of weak transverse perturbations to the dark soliton (A. 1) 

p(x,y,t)=p s (C)-ef(C)e^~ A ^ , 

u(x,y,t) = u s (0-eg(Oe l{m - AT) , 
v(x,y 7 t)=evh((y^r) t 

( = v(x - Kt) , 7] — v 2 y , t — u^t , < e < 1 . 
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-fu' s - gp's + f [« - «»] + Ps (-.</ + ipv 2 h) = -iv 2 Af, 



(A.2a) 



1/2 Ps In J?> I I A £ tt\ i 3f //5g . / 

— 3- (3/ p s + 4/p s ) + j — + 4gu s 



+ ^[p' s (2f-pWf) + 2f'p>:(0 + ff/: 



+4 [/' + g< (u s -«)] + - (pV/' ~ ^/'") - 4i^A S 



ft' = 



(A.2b) 
(A.2c) 



This is an eigenvalue problem parametrized by p and v for the eigenvalue A = A(p; z/) and 
eigenfunction [/, g, h) T . 



Assuming p ^ {0, ±p C utoff } where p C utoff = ^cutoff / v 2 so that the eigenvalue of interest 
is simple, we expand 2 the coefficient functions p s and u s , the parameter n = vl — v 2 , the 
eigenfunction [/, g, h] T and the eigenvalue A in powers of v 2 : 



f = fo + v 2 A + + ■■■ , 

9 = fo + v 2 9i + ^92 + ■ ■ ■ , (A.3) 

h = -ip J (f + v 2 g x + u A g 2 H ) dQ , 

A(p;i/) = A (p) + ^ 2 A 1 (p) + --- . 



Then, (A. 2c) is automatically satisfied to all orders so we only consider eqs. (A. 2a) and (A. 2b) 



2 The two limits, linearization about the soliton and expanding in v, are not interchangeable. 
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which expand, respectively as 



fi-9i + 



f2sech 2 (C) - I) /, 



+ iA fo+p 2 / f d( 



+^ 2 |/2 -32- ^ech 2 (C) (-2 (f[ + g[) + f' Q + 2p 2 J f Q dQ 
-\fo - i/(sech 4 (C)^ + tanh(C)sech 2 (C) (/o - 2 (/i + fli)) 



+iAi/ + iA /i-4/otanh(C)sech 4 (C)+r / // S A 

1 



sech 2 (C) - ^ ) /',, 



8«? 2 - 8/ 2 + 8/ sech 2 (C) tanh(C) (2sech 2 (C) - l) 

+2sech 2 (C) [-4tanh(C) (fg + 2. 9l ) + f»' + Ag[] 
+2sech 2 (C) (6 - 8sech 2 (C)) ft - (2p 2 + l) f Q 

+2f[" - 4g[ + SzAi/o + 8iA Q gi 



(A.4a) 



(A.4b) 



A.l. KP eigenvalue problem. Adding (A. 4a) to (A. 4b) gives 
\ \fd" ~ 4 [(1 - 3sech 2 (C))/ ] ' + 8*A /o + 4p 2 / f d( 



0{» 2 ) 



Differentiating and keeping only leading order terms gives 

Afo = fo" - 4[(1 - 3sech 2 C)/o]" + 8iA /^ + 4p 2 / = . 

This is the KP eigenvalue problem studied in [28]. The unstable portion of the spectrum 
includes one eigenpair 



/o(C;p) = ^2 < ' 



Vl-2p/V3)c 



2-2p/v / 3 + 2yi-2p/v / 3 
V3 



[1 - tanh(C)] 



Ao(p) =-i|y3-2v / 3p, <p < — -^cutoff , 0<^«1 



This eigenvalue is continued onto the positive real line by the eigenpair 

d 2 



MC;p) = ^ <' 



(l-i-v/zp/v^-ijc 



2 - 2p/V3 + 2iJ2p/V3 - 1 



[1 - tanh(C)] 



Ao(p) 



P 



V3 



Pcutoff , < 1/ < 1 



A.2. Perturbed KP eigenvalue problem. Below we determine the correction Ai(p). 
/i is determined in terms of 171 by subtracting (A. 4a) from (A. 4b) to obtain 

V[ - 2ffi - \fo + [sech 2 (C)/o]' +P 2 J fod( = 0{v 2 ) , 
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so that 

fi = 9i + \fo ~ ^ech 2 (C)/ - y J J fod( + 0{v 2 ) (A.5) 
= 9l +f + 0{v 2 ) , / = \fo - ^ech 2 (C)/ - y J JfodC- 

Using (A.5) in eqs. (A. 4a) and (A. 4b), adding the two equations together and differentiating, 
the 0(v 2 ) terms equate to 

£ 9l = - sech 2 (C) [4 fa" - p 2 f + f" + 4/) + # 

+ 2 tanh(C)sech 2 (C) ^ + 3/£" - 4p 2 | f d( + 8/' 

+ 8sech 4 (C) (2% - 4/o + 3/) + V 2 fi - 4iAo/' 

- 8tanh(C)sech 4 (C)/^ + tf(0 - /"" + 40/ sech 6 (C) 
+ 2/"-8iA 1 /5 
= Gf(C;p)-8iAi/^. 

Solvability then determines Ai 

Al(p) = roo . . . , (A.6) 

8J_^/ (C;p)^*(C;pMC 

where h((;p) is the homogeneous solution of the adjoint problem 

= h"" - 4(1 - 3sech 2 <>" + SiA^/i' + 4p 2 /i = . 
Since Aq is either purely real or purely imaginary, the solution of the adjoint problem is 



MC;p) 



/ *(C;p), p<^, 
/o(C;p), 



The integrals in (A.6) can be calculated explicitly 



6V2V% - 3 

For the asymptotic expansion in (A. 3) to be valid, we require Ao(p) ^> ^ 2 Ai(p). This puts a 
restriction on the values of p where the expansion is valid: 

Pcutoff <P<-lr or < p < Pcutoff , < v < 1 . 
Then, the unsealed eigenvalue f2(fc; 1/) in (2.7) has the asymptotic expansion 

Q(k; v) ~ v 3 A (k/v 2 ) + v 5 Ax(k/v 2 ) , 

k < fecutofr(^) , fccutoff(^) < fc = C(i^ 2 ) «1, < v < 1 , 

which is given in (2.13). 

Appendix B. Theorem 2.1. 
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In [30] it was proved that L has exactly one negative eigenvalue which was determined 
explicitly in [8] (L + k 2 utoS /2)f — 0. In addition, it was proved in [57] from general consid- 
erations of linear operators of the form JL where J is skew-symmetric and L is symmetric, 
that the number of eigenvalues of JL with a positive real part is at most the number of neg- 
ative eigenvalues of L. The latter decomposition applies to (2.8), where L = Lq + k 2 /2 is 
symmetric for tgl and J is skew-symmetric. 

Combining these results, for < |fc| < fc C utoff. k E K, L has one negative eigenvalue 
and therefore JL has at most one eigenvalue with a positive real part. By the instability of 
the dark soliton, proven in [30], JL has exactly one eigenvalue with positive real part. There 
is also exactly one eigenvalue with negative real part via the following 

LEMMA B 1 . For fceM, the eigenvalues of JL come in pairs of opposite sign. 

Proof. For any fcel, 

JL = -Ra 3 LJRa 3 . (B.l) 

Let J Li = Tf . Using (B.l) and one of the Pauli matrices (5. 1) gives 

JLRa 3 f = -TRa 3 i . 



It follows that (— r, Ra 3 i ) is also an eigenpair for JL. □ 

On the other hand, for \k\ > fc C utoff, k € M, L has no negative eigenvalues and therefore, 
by Lemma B 1, JL has only purely imaginary eigenvalues. We find numerically and asymp- 
totically in the shallow regime only two discrete, simple eigenvalues for k 6 C\{0, ±fc cu toff }■ 

Appendix C. Proof of Lemma 5.1. 

We make use of the following identity 



JLq = R(T\Lq JRo\ . 
For any k € C, consider (Oq, fo) an eigenpair for eq. (2.7) satisfying 



./ [ L a + ht 2 J +in 
Since Oo is a simple eigenvalue, it follows that 

1 



fn = 



(CI) 



(C2) 



dim kcr 



JL 



-k 2 J 
2 



iQ 



= 1 . 



Therefore, it remains to verify that i?<7ifo ^ s m tne nullspace of [JLq + ^k 2 J + iSlo] \ We 
take the complex conjugate of eq. (C.2) and apply the decomposition (C.l) to obtain 

1 



-Ra 1 L JR(T 1 + -k* J - if2p 







Applying Ro\ yields 



LqJ + —k* J + iQ,Q 



which is precisely the adjoint equation to (C.2). Therefore, we have 



ker 



J[Lq + -k 2 



= span{i?<7ifQ }, k G C. 



By similar arguments with JLq = — o-zLqJoi, one can show that o^fo oc Raxf^ and hence 
spans the kernel of [JLq + \k 2 J + iCloy when k, flo £ R. We use this null eigenfunction 
in our numerical computations whenever k £ (fc cuto ff, oo). 
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